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We calculate the third order nonlinear optical response in the Hubbard model within 
the spin density wave (SDW) mean field ansatz in which the gap is due to onsite Coulomb 
repulsion. We obtain closed-form analytical results in one dimension (ID) and two dimension 
(2D), which show that nonlinear optical response in SDW insulators in 2D is stronger than 
both 3D and ID. We also calculate the two photon absorption (TPA) arising from the stress 
tensor term. We show that in the SDW, the contribution from stress tensor term to the low- 
energy peak corresponding to two photon absorption becomes identically zero if we consider 
the gauge invariant current properly. 

KEYWORDS: Third harmonic generation, Two photon absorption, Gauge invariance, SDW 
insulator 



1. Introduction 

The realization of ultra fast networking through all-optical switching in modern optical 
technology requires advanced optical materials with large third-order nonlinear optical sus- 
ceptibility x^-^H^^^f- !)• Quasi one dimensional (ID) vr-conjugated polymers offer x^'^'^ values 
of 10~^^ to 10~'^ e.s.u. (electronic system of units). The quasi ID Mott insulators such as 
Sr2Cu03, offer x^'^'^ values in the range 10~^ to 10~^ e.s.u.^'^ 

The canonical model for strongly correlated materials is the standard Hubbard model. 
In the spin density wave (SDW) mean field approximation, this Hamiltonian reduces to a 
simple quadratic form that can be diagonalized exactly. Such an SDW Hamiltonian is indeed 
relevant to several groups of organic materials.^' ^ The SDW approximation assumes long 
range order. On the other hand, since such a model is essentially non-interacting, there are no 
vertex corrections in the current loops and the calculation of x^^^ is greatly simplified. This 
simplification allows us to calculate various nonlinear optical processes, including those due 
to the stress tensor. 

Although in one dimensional organic materials, where SDW states have been observed, 
people have already studied the optical conductivity (linear response), nonlinear optical re- 
sponse in this approximation which admits closed form expressions in ID and 2D has not 
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been addressed yet. In this contribution we examine the nonhnear optical response of SDW 
insulators in one, two and three dimensions. Contrary to commonly accepted intuition that 
lower dimensions is equivalent to larger optical responses, we find that among SDW insulators, 
optical response in 2D larger that both one and three dimensions. 

Among ID, 2D and 3D systems with large on-site Coulomb interaction, the ID system 
has the largest optical nonlinearity because of the decoupling of spin and charge degrees of 
freedom.^' In contrast to this, among SDW-ordered systems, the largest third order optical 
response appears in 2D. 

Another purpose of this contribution is to clarify the importance of gauge-invariant treat- 
ment in a simple model. From symmetry point of view, some of optically allowed peaks (such 
as the case of two photon absorption) may become zero, provided there is charge conjugation 
symmetry.^ But, in SDW case, the mean field factorization of the Hubbard model that leads 
to SDW Hamiltonian, breaks this symmetry. However, when dealing with the contribution 
arising from stress tensor terms, we find that the gauge symmetry gives identically zero con- 
tribution to mid-gap peak in TPA. This result implies that the mid-gap peak in TPA is solely 
due to the four-current correlations. This observation is a symmetry property, and as we will 
explicitly show, is independent of dimension. 

This paper is organized as follows: In section 2 we review the SDW mean field treatment 
of the Hubbard model. In section 3 we discuss the choice of gauge and method of calculation. 
In section 4 we discuss the four-current contribution to the third harmonic generation (THG) 
spectrum in ID, 2D and 3D. In section 5 we consider the effects of stress tensor terms that 
come through quadratic couplings of gauge field to the electron system and also through the 
dependence of gauge invariant current to the stress tensor. Finally in section 6, we summarize 
the results. 

2. Model Hamiltonian 

The SDW Hamiltonian is obtained from the Hubbard Hamiltonian in the mean field 
approximation, where the gap is driven by Coulomb repulsion. The Hubbard Hamiltonian is 
written as: 

H = Y1 ^kclcks + U ^(n,-| - n/2)(n,| - n/2), (1) 

ks j 

where cj.^ creates an electron with momentum k and spin s =t, i- The dimension of the lattice 
can be arbitrary, but in this derivation let us focus on ID case. We are interested in half-filled 
case (n = 1). Hereafter we will fix the scale of energy by setting 2to = 1, where to is the 
hoping amplitude. In this paper we also use the system of units in which h = c = e = a = l, 
where a is the lattice constant. The hopping part of this Hamiltonian is characterized by the 
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dispersion = — cos k. The SDW mean field approximation amounts to requiring 

{njs) = ^ + sme''^^, (2) 



where n is the average particle density, and s = it for | and j, respectively, and Q = n. 
Ignoring fluctuations, and dropping additive constant of energy, we obtain the quadratic SDW 
Hamiltonian^^ 

ijSDW ^ Ho-UmY, e'^^ (n,T - n,i), (3) 

j 

where Hq is the tight-binding band part. This can be written in a more compact form as 

^'°^ = EExLw^^x.., (4) 



k s 



where xL = (4.' 4+Qs) = A 4l) and 



^ ,efe + efe+Q ek-€k+Q z tt x 
T-Cka = { 2 — 2 — -Umaa . 

Here A; runs over the half BZ and cr's are Pauli spin matrices. If the perfect nesting property 
ek+Q = — efe holds, we have a much simpler Hamiltonian 

HSdw = _ ,Aa^ (5) 

where the A:— independent gap parameter A = Um is determined by U. 

The unitary transformation ipks = U^Xks, such that UTCkU^ is diagonal, is given by 

..t ( Vk \ I ,2 , I ,2 1 fa\ 

= \ , \uk\ + \vk\ = 1, (6) 

V -4 < ) 

with 

is 

"^= = 71 

—i 



(7) 



£k 



We also need to note the relations: 



2 2 I |2 I |2 ^fe 
'"k-'^k = Ffe - hik = — > 

* * r, 

- = 2ukVk = — , 
^k 

giving 

^k 

U^a'U = +^a' - — (7^ 



(8) 



(9) 
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that imply: 



(10) 



ks 

In ID the density of states for this model is given by (Appendix B) 



(11) 



where wi = \/l + A^. 



One can also consider other quadratic models with gap and coherence factors, such as the 
U = limit of the ionic Hubbard model^^ 



where 'b' stands for band insulator. This Hamiltonian describes simple tight-binding insulator 
in which the gap is due to difference = — = A in site energy, not the Coulomb correlation 
(in SDW A = Urn). This model will differ from SDW insulator in the A replacements 

in their coherence factors. Such a difference will affect the first order responses dramatically, 
but it is easy to see that third order optical response of this model is identical to SDW model. 

3. Choice of gauge 

The coupling of electromagnetic field to matter can be described in two gauges. One is 
a gauge in which vector potential is zero, but the scalar potential is non-zero and given by 
Aq = — E.r. In this gauge, the electric field of radiation couples to the dipolc moment of 
electrons, and hence one needs the matrix elements of the position operator r to calculate 
the response of the matter to electromagnetic perturbation. Working in this gauge is suitable 
for molecules and small clusters. Because of the r operator, the calculations in this gauge are 
sensitive to boundary conditions. Moreover, for periodic boundary conditions, one has problem 
in choosing the origin of the r coordinate. Therefore this gauge is ill-defined in thermodynamic 
limit, and sensitive to boundary conditions. 

On the other hand, we have an alternative choice of working in a gauge in which scalar 
potential is zero, while the vector potential A is non zero. In this gauge the coupling between 
the external field and electrons at the first order is via the current operator: j.A. In second 
order the gauge field couples to electrons via the stress tensor operator as A.r.A, etc. Working 
in this gauge is actually equivalent to Peierls substitution, and hence we call it the Peierls 
gauge. 

Without taking into account the effect of nonlinear couplings, and nonlinear dependence 
of the gauge-invariant current on the vector potential, the response function at n'th order is 
given by: 



= -to J2 {4be±i + bl^T^ae) + Yl 




(12) 




X 



(")(J^;a;i,..., 



ne- 




Xj" 



) 
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where n— current correlation function is given by 




(14) 



Here Tc is the time-ordering operator along the Keldysh path^^ and j(r, t) is particle current 
operator. 

In general a m— point Keldysh Green's function has a tensor structure due to two time 
branches. A unitary transformation to the "physical" representation will give Gap-y... where 



Here a stands for "advanced", while r means "retarded". The optical experiments measure 
the fully retarded components which are given by nested commutators. On the other hand, 
the sum of nested commutators is generated by Garr...,G'tc., where only one of the indices 
is equal to a and the rest are equal to r index. ^ As will be shown in appendix A, due to 
the commutators and appropriate 9 functions, this component of Keldysh Green's function is 
very special, in the sense that we do not really need to get into Keldysh machinery in order 
to calculate the nonlinear optical response. As is shown in appendix A, the fully retarded 
Garr...r Component can always be calculated within the framework of equilibrium quantum 
field theory. We first calculate the time ordered expectation values and then analytically 
continue the result to ensure the correct behavior of the poles. If we need other components 
of Keldysh Green's functions (the fluctuation functions) that involve anti-commutators, and 
are related to noise spectroscopy, Keldysh formulation becomes inevitable. 

Therefore, in the case of optical response, we can forget the two time branches, and denote 
the response function at say, third order by (jjjj), keeping in mind that this is an ordinary 
time ordered expectation value and should be analytically continued according to prescription 
of appendix A. 

Another important point is that the four-current scheme for calculation of the optical 
response is reliable far from zero frequency. Therefore, we do not have to worry about the 
so called zero frequency divergence (ZFD) in our calculations.^^ To remove the unphysical 
ZFD one has to calculate a few more correlation functions, but as far as the behavior near 
resonance region is concerned, the (jjjj) is sufficient. To appreciate this point, let us look at 
the first order response in two gauges (// is the dipole operator and is proportional to r; j is 
current operator, and proportional to r): 



a,/3,7 G {a,r} 




(15) 




(16) 
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where a factor of 2tt5{uj + uJa), with tUu = cui does also multiply the right hand side. We 
integrate by parts with respect to both time variables r and ri to obtain: 



Jo Jo 

Jo 



lUJl 



^dre'^^ [(/x(r)M/3)) - (MtHO))] 



— - f^dn e-i-i MPWi)) - (/x(0)j(ti))] 

lOJl ILJ Jq 



ILOl lUJ Jq 

+ ^^^^ 

where we have used the fact that both ui and uj are bosonic Matsubara frequency, e^'^'^^ = 1. 
Hence we see that the relation x^^(a;;a;i) = ■^^Xjji'^'i^i)^ often used to relate the response 
in two gauges is not quite correct. Inclusion of the boundary terms (first two lines of the 
above equation) compensates the frequency denominators appearing in [jj) scheme. Similar 
consideration applies to higher order correlations such as {jjjj)- 

If we keep nonlinear coupling between the gauge field and electrons that comes through 
the stress tensor operator r, and also dependence of gauge invariant current to powers of 
gauge field A, a simple perturbation theory gives expressions of the type {jjt) that again 
appear in fully retarded combination of nested commutators. The simplicity of SDW and our 
band model allows us to consider these terms as well. 

4. Four-current response 

The particle current operator is given by 

J = ito ^ (4+1'=^ ~ 4^i+i) (18) 
e 

The momentum space representation of this operator for the band and SWD insulators is 
given by: 

J = XI xlslk(^^Xks, Ik = 2*0 sin k (19) 

ks 

which in terms of new fermions is 

/^^ = E^Lf^-^-^-^V^^ (20) 

ks \ £k / 

The coefficient of in the above expression describes mira-band transitions, while the coef- 
ficient of causes mier-band transitions. Let us define the corresponding coefficients by 

g£k = sink cos k, hsk = sAsin/s, (21) 

The THG susceptibility corresponding to photon frequency v (or iu in imaginary time) is 
given by 



6/19 



J. Phys. Soc. Jpn. Full Paper 

A few comments are in order: The trace operator arises from closing the current loop, and 
enormously simplifies the calculations. This trace is taken with respect to the indices of Pauli 
matrix. But since Pauli matrices along with unit matrix have a closed algebra, and that only 
the unit operator survives the trace, the sums of 16 terms in the case of first order response, 
and of 256 terms in the case of third order response are considerably simplified. 

We need one more physical considerations to further simplify the calculations. At half 
filled situation of interest to us, the intra-band terms at zero temperature do not contribute to 
optical absorption. In calculating general expectation value of say (ABCD), where A, B, C, D 
can be any operators contributing to the coupling of light with matter, we are interested in 
dominant processes in which terms containing g (intra-band matrix element) in the rightmost 
and leftmost operators A and D do not contribute. Similarly in B and C operators the h term 
should be dropped. 

To calculate the time ordered product within the framework of equilibrium quantum field 
theory, we use the Matsubara technique. The summation over fermion loop frequency iun can 
be done with standard contour integration techniques. Also, if we ignore the momentum from 
the incident light, the momentum k running in the current loop must be the same for all 
fermion propagators. After contracting various fermion spinors to get the appropriate Greens' 
functions we obtain 

^OjJj)=Eix (22) 

8(44 + ;/')A2 cos(fc)^ sin(fc)4 

(2£fc - 3i^)(2cfc + 3i/)(2£fc - 2y){2ek + 2p){2ek - v){2ek + v) ' 

The analytic continuation v ^ v + iO+ is implicitly understood. 

It is also important to note that, as far as the behavior of response functions near the 

resonance is concerned (that is away from = 0), four-current response functions, (jjjj) give 

the same qualitative features as that of four-dipole response functions, (///n/z/n), where // is 

the dipole moment operator. To appreciate this point, let us go back to the expression of 

the current operator for the SDW insulator. The inter-band matrix element for the current 

operator is jk = sAsinA;/£fe, that via the equation of motion for the dipole moment operator 

jjL, would imply that the dipole matrix elements are: 

/Xfc = sAsinA;/2e^. (23) 

Similar procedure that lead to equation (23) gives: 

= ^^y. (24) 

k ^ 

(4£^. + t/2)A^ cos(fc)^ sin(fc)-^ 

(2efe - 3z/)(2efe + ■iv){2ek ~ 2v){2ek + 2v){2ek - p){2ek + v) ' 
Now we can see that up to a numerical factor of the order of unity, the "near resonance" 
behavior of both expression is the same. Because after analytic continuation, the imaginary 
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Fig. 1. Plot of imaginary part of x™*^(t^) vs. v in the SDW model in = 1,2,3 dimensions. The 
unit of energy is 2io ~ 1, and in all calculations we have assumed A = 0.3. The gap is given by 
Eg — 2A ~ 0.6. We see that the response in D~2 is stronger than D=l,3. Values of D=l,3 are 
magnified by factors of 10 and 100, respectively, for eye assistance. 



parts give delta functions peaking near ~ e^. Therefore i^'^ef. in the denominator of (jjjj) is 
of the same order of magnitude as in {^fi^jj) expression. Note that four-current formula 
makes sense near the resonance conditions only. It gives the unphysical divergence in the 
static limit ^ 0. Fixing this problem as discussed in equation (17), requires calculation of a 
few more expectation values, but we are not interested in this limit here. Since we are dealing 
with a gapped situation in which the frequencies of interest are of the order of gap and zero 
frequency is avoided. 

The imaginary part of THG susceptibility can now be written as 



27, f3u 



EvrA^ cos^ k sin'^ k ^ . 
—6 ^(^k 

k 



(25) 



where D is the dimension of space. In ID this integral can be evaluated as follows: 



/i(^) 



dep{e) 



^A2(e2_A2)(u;2-e' 



2\2 



24|i/|5 

Using the Kramers-Kronig relation 



12e6 



-5{e 



A 



(26) 



(27) 



h' — U) 

one can obtain a closed form expression for the real part of the above retarded susceptibility 
in terms of Elliptic functions of various kind. 

In 2D still we can obtain closed form results bye transforming to tight binding coordinates 
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(appendix B) 

+ 9^1^ (^) f^^^'^' " ^^^^ + • ^^^^ 
In this expression = u'^/A — A^. Here E and K are elliptic functions of second and first 

kinds, respectively.^^ The second term near the band edge behaves like — Alog(z//2— A) 

which remains finite, while the first term give a divergent contribution of the form /2(z^/2) ~ 
1/1/1^/2 - A. The main THG peak is due to /2(3i//2) ~ l/-^3i//2 - A which occurs at 
1/ = 2A/3, with 2A being the excitation gap. 

In three dimensions fsii') and hence can be calculated numerically. Figure 1 

shows Qx™'^{^) for -D = 1, 2, 3 and the gap parameter A = 0.3. The 3D result is magnified 
100 times (dotted line), and the ID result (dashed line) is magnified 10 times. In the SDW 
system, inverse square root divergence in 2D has no counterpart in ID and 3D. To understand 
this, notice that, one can always replace a k integration with a energy integration weighted 
by DOS, and so such an enhancement in 2D compared to ID and 3D can be traced back to 
the nature of singularities of DOS (appendix B). This figure demonstrates that if we have a 
insulator in which gap is due to SDW type of order, the phase space effect (DOS) along with 
SDW coherence factors (u^, Vk) generate larger nonlinear response in 2D compared to ID and 
3D. 

5. Stress tensor terms 

The simplicity of SDW model, along with possibility of obtaining closed form results in 
one and two dimensions allows us to investigate the role of stress tensor terms in nonlinear 
optical processes. The first place that stress tensor appears is via the coupling of external field 
to matter at second order which is r.A.r. Then the gauge invariant current (particle current 
plus a correction from external field) involves the stress tensor itself, and to first order is given 
by:i5 

jm^jm_ jranj^^^ (29) 

Therefore the third order response of gauge invariant current involves the response of 
both 'particle- current operator {f^ above) and also the response of stress tensor operator 
(^mn above). The third order response due to particle-current is 

(jrj) + (rjj). (30) 

while the third order response due to r operator is 

W (31) 
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These expressions are of the type given in (A- 10). The first two corresponds to A,Q ^ j, B ^ 
T, and A ^ T, B,Q ^ j, respectively, while the third one corresponds to A, B ^ j, Q ^ t. 
To obtain the appropriate matrix elements, let us begin by writing the stress tensor in ID 

as 

^ = -ioJ2 {4isCn+ls + 4i+lsCns) , (32) 

ns 

or equivalently r = Ylks xls^k'^^Xksj that eventually becomes 

^SDW Y^,/,t f^k^z /, /qq\ 

^ =2^V'fe.^-^ JV'fcs. (33) 

The matrix elements we need are given in table I, where s is ±1 for | and J, spins. 





jinter^intra^inter 


jinter^intra^inter 


SDW 


sin^ /ucos^ k/el 


sA^ sin^ /ucos^ k/el 



Table I. Inter band and intra band matrix elements in the SDW model. 



respectively. Since the matrix elements involve both sin and cos, in ID the transitions at zone 
center and zone boundary are suppressed, and hence there is no divergence. But in 2D again 
SDW response will be divergent. 

Let us look at the two photon absorption (TPA) contribution arising from r operator 
carefully. In ordinary case of (jjjj), the TPA corresponds to uji = L02 = —"^3 = J^- In the case 
of stress tensor TPA corresponds to uji = UJ2 = i^- If we denote the denominators of the first, 
second and third lines of the expression (A-10) by ^1,^2, 4; respectively, after decomposing to 
partial fractions we have 

^^'^ 2£k{v + i-q - Ek) 2efc(i/ + ir/-2£fc)' ^^^^ 

= 2ek{v + ir^-2eky ^^^^ 
where & {y + ir] ^ —v — ij]) term also is present in the time ordered correlation function. 
Putting in the matrix elements, the TPA susceptibility becomes 

(y"**^'")V"^''-^M4 + ^2-4] (36) 

+ y-t^v-f-v'"*^-- [(4 + 4) + (4 + ^1) - {h + h)\ , 

where we have not simplified the second line deliberately to emphasize the role of gauge 
invariant current. Obviously what we measure is the gauge invariant current. The minus sign 
in third terms of the each line in the above equation comes from the minus sign in equation 
(29), which is essential for the cancellation that takes place in the second line of the above 
equation, leaving a contribution proportional to oc ^3. 
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Now let us concentrate on the first line of this equation. Here the matrix element have 

opposite signs for up and down spins (which ultimately comes from the coherence factors of 

the SDW state). Therefore their contributions identically cancel each other, and we are left 

only with a term coming from denominator of £3 only, which peaks at the gap 2 A. In ID as 

we mentioned, the presence of sin and cos suppresses the transitions at band edge and zone 

center. So let us calculate this term in 2D in the SDW insulator: 

// \\SDW ttA^ ■t-^ siia^kxCos^kx 
{{jxJxTxx)) = -7- 2^ 5{v/2-ek) 4 

"■X ■t"'y 

ttA^ f f dxdw sin^.7:cos^x ,„ , 
= -—JJ^^^'^''^'-'^ (37) 

A"^ f f dA(i£ sin^ xcos^ X , 
= -^1]— ? ^^""'^-'^ 



-4A^^^- ^"^^ - 16- 



where in the last line again delta function picks up the value of A = — A^, P = 

1 — as in the appendix B. The logarithmic divergence of the first term gets suppressed 
by A ~ -x/v — 2 A factor, but the second term still shows inverse square root divergence in 2D. 
Therefore in the case of stress tensor terms too, the 2D SDW systems offer larger response 
than ID. 

We would also like to emphasize that, the requirement of gauge invar iance cancels the 
singularity at A = Eg/2, in r's contribution and leaves us with a square root singularity in 
2A = Eg in 2D SDW systems. In other words, in SDW insulator the peak at Eg/2 may solely 
be due to four-current correlations and the contribution from stress tensor to this peak is 
identically zero. 

6. Summary and conclusions 

In conclusion, we have calculated the nonlinear optical responses in the SDW insulator in 
which the gap is due to on-site Coulomb repulsion. The linear response of SDW model has 
the characteristic inverse square root divergence in ID which due to larger amount of nesting 
is further enhanced by a logarithmic factor in 2D and is entirely suppressed in 3D. The THG 
spectrum of SDW model has no divergence in ID, but diverges as inverse square root in 2D 
and becomes finite again in 3D. Therefore the optical responses (linear and nonlinear) of SDW 
insulators is maximal in 2D as a function of dimensionality. 

The model calculations presented in this work suggests that nesting as a possible mech- 
anism of nonlinearity enhancement which works best in 2D, contrary to common intuition 
that lower spatial dimensions are better for nonlinear optical materials. This mechanism does 
not work in ID. It was found that in ID such a enhancement can arise from the spin-charge 
separation.^' 

The simplicity of the quadratic model treated in this investigation allowed us to exactly 
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calculate the contributions of stress tensor in nonlinear optical response. We showed that (i) in 
TPA measurements in the SDW systems, the structure at the mid-gap {Eg/2) has essentially 
no contribution from the stress tensor term, and (ii) when these contributions are non-zero, 
have comparable effect to that of usually considered four-current terms. Stress term has even 
parity and in nonlinear processes can lead to dipolc-forbidden transition. This hints to the 
importance of gauge invariant treatment of currents in nonlinear optics, which is usually 
neglected in the literature. 
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Appendix A: Analytic continuation 

In nonlinear response theory, wc need to calculate fully retarded expectation value of 
nested commutators. For example in the case of nonlinear dielectric response of an electron gas 
to a fast moving ion, the nested commutator of density operators at different times appears. 
In the general theory of nonlinear response, wc might have the fully retarded combination 
of nested commutators of arbitrary operators. These operators arc determined by nature of 
coupling (linear, quadratic, etc.) between the system and the external perturbation, and also 
the observable being studied. For example the second order coupling of the electromagnetic 
field to matter via the stress tensor operator r, leads to a fully retarded current response of 
the form ([j, [j, r]]). The general framework to study this type of expectation values is the 
non-equilibrium quantum field theory. However, within the standard formulation of quantum 
field theory, with appropriate analytic continuation, one can obtain these kind of expectation 
values that correspond to Garr...r component in Keldysh Green's function language. In optical 
measurements always this kind of expectation values appear. In noise spectroscopies the other 
components of Keldysh Green's function appear that can not be treated as straightforward 
as Garr...r Component, and use of Keldish formulation becomes necessary. 

For the problem in which external perturbation couples linearly to the system, through 
operator Aj , and quadratically through B^i , where we are interested in variations of quantity 
Qi, one can write the retarded response at third order as: 

^^M^;h,t2) = +eit - h)9{h - h){[A,{h), [Bkiih), Q^m) 

+ 9{t - h)0{h - t2){[Bki{t2), [Aj{h), Qiit)]]} 

+ {j^k^l^j) + {j^k^l^ jf. (A-l) 
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On the other hand, the time-ordered product of these three operators is given by: 

cl>fjkiit;h,t2) = {TAj{ti)Bki{t2)Qi{t)) = 

+ 0{tl - t2)e{t2 - t){Aj{ti)Bkl{t2)Qi{t)) 

+ eit - t2)e{t2 - ti){Qi{t)Bkiit2)A,ih)) 
+ e{t2 - h)d{h - t){Bki{t2)Aj{ti)Qi{t)) 

(A-2) 

+ 9{t - ti)9{ti - t2){Qi{t)Aj{h)Bki{t2)) 

+ e{ti - t)e{t - t2){Aj{ti)Qi{t)Bki{t2)) 

+ e{t2 - t)e{t - ti){Bki{t2)Qiit)Aj{ti)) 

+ permutations. 

The signs are all positive, since operators A, B, Q are quadratic in fermion operators. 

Now let us expand the commutators in definition of retarded expectation value to obtain 

(pfjkii't;ti,t2) = 

+ eit - t2)e{t2 - ti){Aj{ti)Bki{t2)Qi{t)) 

- e{t - t2)e{t2 - ti){Aj{ti)Qi{t)Bki{t2)) 

- e{t - t2)e{t2 - tl){Bkl{t2)Q^{t)A,{t^)) 

+ e{t - t2)e{t2 - ti){Qi{t)Bki{t2)Aj{ti)) 
+ eit - ti)eiti - t2){Bki{t2)Aj{ti)Qi{t)) 

- 9{t - ti)9{ti - t2){Bki{t2)Q^{f:)A,{ti)) 

- 9{t - ti)9{ti - t2){Aj{ti)Qiit)Bki{t2)) 
+ 9it - ti)9{ti - t2){Qiit)Aj{ti)Bkiit2)). 

Now using the identity 

0{t - t2)0{t2 - h) + 0{t - ti)0{ti - t2) = 9{t - t2)9{t - h), (A-3) 
the sum of second and seventh terms simplifies to 

-9{t-t2)9{t-ti){AQB), (A-4) 
while the sum of third and sixth terms becomes 

-9{t-t2)9{t-ti){BQA), (A-5) 
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so that we obtain 

+ e{t - t2)e{t2 - h){Aj{ti)Bki{t2)Qiit)) 
+ e{t - t2)e{t2 - h){Qi{t)Bki{t2)Aj{h)) 

+ e{t - h)e{h - t2){Bki{t2)Aj{h)Qi{t)) (A-6) 

+ e{t - ti)e{ti - t2){Qi{t)Aj{ti)Bki{t2)) 

- e{t - h)e{t - t2){Aj{ti)Qi{t)Bki{t2)) 

- e{t - h)e{t - t2){Bki{t2)Qi{t)Aj{h)). 

Now apart from the 6 function that determines the analytical structure, the time ordered and 
retarded nested commutators have similar structures. Therefore one can obtain the expecta- 
tion value of fully retarded nested commutator by appropriate analytic continuation of the 
corresponding time ordered one. So, let us obtain the spectral representations of the retarded 
and time ordered expectation values and compare them. 
Our conventions for Fourier transforms are 

/ + 00 J I- + 00 

^e-'^'f{iv) ^ /» = / dte+'^*f{t), 
-oo J — oc 

diO f duJl I d,LO\ ,t . t , t 



X(c.;cci,c.2) =J^J^J ^e-V-^^^e-='^x(i;ii,i2). (A-?) 

Note that the operator -6(^2) acting at time ^2 couples the second power of the external 
perturbation to the system. Therefore, in principle it could involve two frequencies uj2,uj3. But 
since the external fields are supposed to act at the same time, we always have the combination 
L02 + UJ3- Therefore here we have dropped the ^3 in our calculations. Using the representation 

— - , (A-8) 
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for the step function, after some algebra we get 

a.b 



while 



+ 



+ 







+ Em - ir]){'^i - 


Eoa 


- ifl) 




+ EQa + iri) {uji - 
BllAfQf 




+ irf) 




+ Ebo - if]) {uj2 - 
Q^AfBfi 




-irf) 




+ EQa +ir]){^2 - 
AfQfBfi 


EbQ 


+ itf) 




- Em + ir]){uJi - 
BOfQfAf 


Eoa 


-ir]) 



{uj2 - Eoa - i'n){ux - Ebo + iri) ' 



a,b 



+ 



+ 



AfBl\Qf 




+ Ebo + iri){ui - 

QOaBabAf 


Eoa + iri) 




+ Eoa + iri){ijJi - 
BOfAfQf 


Ebo + iv) 




+ Ebo + iri){L02 - 


Eoa + iv) 




+ EQa + iri)(u2 - 
AfQfBfi 


Ebo + iv) 


{0J2 


- Ebo + iri){ui - 
Bl-QfAf 


Eoa + iri) 



(A-9) 



(A- 10) 



{(ji)2 - Eoa + iri){ui - Ebo + ir)) ' 
where we have defined Wo- = uj2 + oji- We see the exact parallelism between the time ordered 
and fully retarded expectation values. The important difference is due to the nature of step 
functions that give rise to Un + irf (n = 1,2, a) structure in the retarded function. This 
is what was expected from causality imposed by appropriate 6 functions in fully retarded 
one. However, as a consequence of having 6 functions along with commutators, equal time 
contractions in diagrammatic perturbation theory do not contribute, and should be excluded. 

We can also write the above result in a more compact form if we note that the frequencies 
are associated with operators as {A, loi), {B, L02), {Q, —uJa)- If we ignore the iij factor, the above 
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result can be written as 

j^Oa j^ab ^bO 



where V stands for all different permutations of (^4, cui), {B, UJ2), (Q, —i^a)- At the end we must 
remember to use the appropriate h] factors to ensure LOn + ir] structure. One can easily see that 
this prescription gives the correct spectral representation in case of two-current correlation 
with A = B = j: 

We can derive a similar prescription for higher order correlation functions of fully retarded 
nature in a straightforward way: 

a^b,c 

(A-12) 

j^Oa^abfjbcQcO 

Appendix B: Tight binding coordinates 

In this appendix, we denote the and ky coordinates in the reciprocal space by .x.y for 
convenience. Since constant energy surfaces cos cosy appear very frequently in calcTilations 
related to 2D tight binding systems, it is useful to define a natural orthogonal transformation, 
(x,?/) ('^lO' that constant coordinate surfaces correspond to constant energy surfaces. 
The first coordinate obviously must be 

A = cosx + cosy. (BT) 

In order to guess an appropriate form for the second coordinate ^, we require the constant ^ 
surfaces to be orthogonal to the constant A surfaces, that is 

oc VA = — sin a; e.^ — siny Cj,, (B-2) 

where is a 'velocity' tangent to the constant ^ surface. This equation implies that 

dx . dy . X 
— = —smx, — = — smy, (B-3) 
dt ' dt ^ ^ 

the division of which gives 

dx dy 



din (tan(x/2)) = din (tan(y/2)) . (B-4) 
Sin X sm y 

Integrating the above equation gives tan(y/2) = const x tan(x/2), where we define this con- 
stant to be tan^: 

tan^ = tan(y/2) cot(x/2). (B-5) 

Equations (B-1) and (B-5) imply that 

dxdy = JdXd^, J = . ^ , (B-6) 

v/l-/3cos2(2^) 
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with (3 = 1- AV4. As a cross check for this formula, one can calculate the area of the 
BZ (J dxdy = 47r^) in the new coordinate system. A straightforward numerical integration 
reassures us that the above Jacobian is correct and the intervals < ^ < 27r, — 2 < A < 2 
count the original BZ only once. 

Also the inverse transformation is given by 

A J-^ -I 

^°^"=2+^^' ^^-'^ 

where the Jacobian J is already defined in equation (B-6). The variable ^ can be interpreted 
as an angle. In fact near the F point where x, y and hence, tanx ~ x, tany ~ y, we 
can actually see that tan^ ~ y/x and hence near to F point ^ can be identified with the 
polar coordinate 0. In this limit also we can write A 2 — (x^ + 2/^)/2 = 2 — r^/2. Therefore 
our coordinate system is a natural extension of polar coordinates (r, 0). Near the zone center 
the energy contours become circular, but near the Fermi energy A = 0, the contours are 
rectangular. 

As an example of the application of this coordinate system, let us calculate the exact DOS 
for SDW systems in 2D. All we have to do is to switch to the new coordinate system so that 



The A integration can be performed by changing to the new variable e = VA^ + A^ that 
gives dX = ede/X. Hence the final result is 




where K is the elliptic integral of first kind,^^ and the restriction < A < 2 translates to 
< u < W2-, with W2 = 4 + defining the upper band edge. This result in the limit of 
A — >^ reduces to the appropriate result for the tight binding bands. In this limit we have a 
\og{v) singularity at the Fermi surface. But for A 7^ 0, the nature of singularity near the lower 



band edge for 2D SDW system is log(z/ — /S) / \/ v — /S.. In ID, the logarithmic contribution is 
absent. Hence the lower band edge in 2D spin density wave systems offers more DOS than ID 
case. We see in the text that this simple observation has profound implications on third order 
nonlinear ity in SDW systems. 

Calculation of DOS in ID is much easier. We need to calculate the following: 

pH= Y.^{ek-u) = j depo{e)5{e-u). (B-U) 

k 

The DOS corresponding to gapless situation A = 0, is given by po{u) = — w^. Now we 

use the relation e = — A^ between energies e and e corresponding to A = and A 7^ 
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situations, respectively. This change of variables gives: 

f ede 5{e — uj) \uj\ 1 fB 12) 

where wi = Vl + A^. For SDW insulator, the gap parameter A = —Um is determined by 
Coulomb correlation. 
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